Recent advances in statistical software1 have enabled public health researchers to fit multilevel models to a variety of outcome variables. Multilevel models facilitate inferences regarding unexplained variability among randomly sampled clusters of units (e.g., hospitals) in outcomes of interest and identify covariates that explain the variance in a given outcome at each level of a particular data hierarchy (e.g., patients within hospitals).2,3 Models with random intercepts enable researchers to accommodate correlations within higher-level units resulting from longitudinal or clustered study designs, and models with random coefficients enable researchers to identify higher-level covariates that explain between-cluster variance in relationships of interest.2,3
Public-use survey data sets collected from large national samples, such as the National Health and Nutrition Examination Survey, also have become widely available.4 The samples underlying these data sets are often “complex” in nature for 2 reasons: (1) the use of stratified multistage cluster sampling to increase sampling and cost efficiency and (2) unequal probabilities of selection from target populations for sampled elements, often as a result of oversampling of key subgroups (leading to the need to use weights for generating unbiased population estimates). Secondary analysts can accommodate these design complexities statistically by using “design-based” analyses, which ensure that population inferences are unbiased with respect to the sample design.4 However, these design-based approaches generally do not enable the types of cluster-specific inferences afforded by multilevel models,2,3 and researchers are now considering multilevel models for complex sample survey data.
Multilevel modeling represents a “model-based” approach to survey data analysis, in which dependencies in the data introduced by complex sampling features are generally accounted for by sound specification of the underlying probability model.5,6 Advocates of this approach argue that any information contained in the sample design features should be accounted for in the model specification, making the sampling uninformative.5 However, analysts may not have access to covariates capturing all of this information. In this case, the use of weighted estimation when fitting multilevel models provides some protection against potential biases introduced by informative sampling.6 Informed by recent methodological and computational developments in this area,1–3,6,7 we show that changes in inferences are possible when fitting multilevel models to complex sample survey data and ignoring the sampling weights.
We analyzed data from the 2013 Medical Monitoring Project HIV Provider Survey, sponsored by the Centers for Disease Control and Prevention, for which a probability sample of HIV care providers was selected from outpatient HIV care facilities in 16 states and Puerto Rico.8,9 Briefly, the provider survey followed a 2-stage probability-proportionate-to-size sample design, first sampling states and territories and then HIV facilities and selecting all providers within a facility. Unbiased estimation of multilevel model parameters requires the use of weights at all levels of a given data hierarchy,7 so we used previously calculated sampling weights adjusted for nonresponse at the facility level and inverses of estimated response probabilities at the provider level.
We focus on only facilities with multiple responding providers and include covariates that are both theoretically relevant for the dependent variables described later in this article and related to the sampling weights (e.g., an indicator of the provider serving more than 200 patients). Details about computation of the Medical Monitoring Project sampling weights for both providers and facilities are available on request.10 We scaled the final provider-level weights to sum to the sample sizes within each facility. A failure to do this would overstate actual sample sizes within each higher-level unit (facility), possibly resulting in biased estimates of model parameters.2,3,7
We fit multilevel logistic regression models to 2 binary dependent variables, indicating whether the responding provider delivered adequate drug use risk reduction and sexual risk reduction services to patients (defined as delivering approximately 70% of recommended risk reduction services to most or all of the patients). The models included random intercepts to capture between-facility variation in each proportion, in addition to fixed effects of several provider- and facility-level covariates of interest. We fit these models with the new GLIMMIX command11 in SAS/STAT version 13.1 (SAS Institute, Cary, NC), which can fit multilevel models to complex sample survey data. Identical results can be obtained with the new svy: melogit command in Stata version 14 (StataCorp LP, College Station, TX).
We did not test whether the parameter differences in the weighted and unweighted models were significant,12 but we did observe several shifts in inference when using weighted estimation (Table A; available as a supplement to the online version of this article at http://www.ajph.org). In both models, the intercept became more negative and significant, suggesting that the probability of using adequate risk reduction was being overstated for the type of provider represented by zeroes on all of the covariates (which may not be entirely meaningful in all models). For drug risk reduction, the coefficient for delivering care in a language other than English became nonsignificant. For the sexual risk reduction outcome, the male provider coefficient became significant, and the Black provider, nurse practitioner, and integrated team effects became even stronger. Finally, the estimated variability of the random facility intercepts was clearly being overstated when ignoring the weights, and the weighted models explained more of the variance in the outcomes at each level.
The weights at each level were clearly informative about the parameters defining these models, and ignoring them in analysis would have led to erroneous inferences with respect to the sample design used. Notably, these results held despite the inclusion of available covariates related to the sampling weights in the models. In practice, covariates used to compute the weights or the weights at each level of the data hierarchy may not be available to the public, making appropriate design-adjusted estimation of multilevel models difficult or impossible. We encourage analysts fitting multilevel models to survey data to carefully examine the variables available for weighted estimation in these data sets, make use of the powerful software1–3,11 that has been developed in this area, and (when possible) examine whether weighted estimation or adjustment for covariates related to the weights affects their inferences.